Simulations of thermal Bose fields in the classical limit 
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We demonstrate that the time- dependent projected Gross-Pitaevskii equation (GPE) derived 
earlier [Davis et al, J. Phys. B 34, 4487 (2001)] can represent the highly occupied modes of a 
homogeneous, partially-condensed Bose gas. Contrary to the often held belief that the GPE is valid 
only at zero temperature, we find that this equation will evolve randomised initial wave functions 
to a state describing thermal equilibrium. 

In the case of small interaction strengths or low temperatures, our numerical results can be com- 
pared to the predictions of Bogoliubov theory and its perturbative extensions. This demonstrates 
the validity of the GPE in these limits and allows us to assign a temperature to the simulations 
unambiguously. 

However, the GPE method is non-perturbative, and we believe it can be used to describe the 
thermal properties of a Bose gas even when Bogoliubov theory fails. We suggest a different technique 
to measure the temperature of our simulations in these circumstances. Using this approach we 
determine the dependence of the condensate fraction and specific heat on temperature for several 
interaction strengths, and observe the appearance of vortex networks. Interesting behaviour near 
the critical point is observed and discussed. 

PACS numbers: 03.75.Fi, 05.30.Jp, ll.lO.Wx 



I. INTRODUCTION 

The observation of Bose-Einstein condensation (BEC) 
in dilute alkali gases ||, ||, ^ heralds a new era in the 
study of quantum fields. It offers a unique opportunity 
to carry out experiments in the laboratory for which the- 
oretical calculations beginning from a microscopic model 
of the system are tractable. However, such calculations 
are fraught with difficulties at finite temperatures. While 
equilibrium perturbation theories have had much success 
[g, ^ 1^ dynamical calculations often require severe ap- 
proximations to be made. 

In Ref. we developed an approximate formalism 
to describe the dynamics of a thermal Bose condensate 
based on the Gross-Pitaevskii equation (GPE). This de- 
scription is valid when the low-lying modes of the sys- 
tem are classical, satisfying the criterion Nk ^ 1. This 
is analogous to the situation in laser physics, where the 
highly occupied laser modes can be well described by 
classical equations. We proceeded by dividing the field 
operator into a classical region represented by a wave 
function ip{'x.) describing the condensate and its coher- 
ent excitations, with the remainder of the field described 
by the quantum operator ?7(x). We derived an equation 
of motion for ')p{'x.) that we called the finite temperature 
Gross-Pitaevskii equation (FTGPE). 

The FTGPE is a rather complicated equation, how- 
ever, and in Ref. ||] we briefly described the first re- 
sults from the simpler projected Gross-Pitaevkskii equa- 



tion (PGPE) obtained by neglecting the operator fi{x.). 
These results demonstrate that the GPE alone can rep- 
resent thermal Bose gases. In this paper we elaborate on 
these results and describe our method in more detail. We 
also consider the effect of strong particle interactions on 
the thermal distributions and investigate the appearance 
of vortices in our simulations. 

The use of the dynamical GPE at finite temperature 
was originally proposed by Svistunov, Kagan, and co- 
workers |l^, |l^. Despite this suggestion first 
appearing in 1991, there have been relatively few numer- 
ical studies based on this approach. Damle et al. have 
performed calculations of the approach to equilibrium of 
a near ideal superfluid jl4j , while Marshall et al. ||T^ car- 
ried out a qualitative study of evaporative cooling using 
a 2D GPE. References |l|. 



17 



m m, pi also use 



'Electronic address: mdavis@physics.uq.edu.au 



References 

classical methods to represent thermal Bose-condensed 
systems. Similar approximations to other quantum field 
equations have been successful elsewhere [ p2[ . 

This paper is organised as follows. In Sec. II we give a 
brief derivation of the finite temperature Gross-Pitaevskii 
equation. In Sec. Ill we describe and justify the simpli- 
fication of the FTGPE to the projected Gross-Pitaevskii 
equation, before describing the simulations we have car- 
ried out in Sec. IV. Section ^ presents the qualitative 
evidence that the simulations have reached equilibrium, 
while Sec. VI carries out a quantitative analysis of our 
numerical data. Section VII discusses the behaviour of 
the condensate fraction, specific heat, and vorticity of 
the system with temperature, before we conclude in Sec- 
tion IVIIll. 
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II. OUTLINE OF FORMALISM 

A full derivation of the FTGPE and a discussion of the 
physics described by each of the terms can be found in 
Ref. §. Here we outline the derivation beginning with 
the equation of motion for the Bose field operator 

= Hsp^i^) + C/o4'^(x)*(x)*(x), (1) 

where Uq = ATrh^a/m is the effective interaction strength 
at low momenta, a is the s-wave scattering length, and 
m is the particle mass. H^p is the single-particle hamil- 
tonian defined by 



2m 



+ Vt,ap(x), 



(2) 



where Vtrap(x) is the external trapping potential, if any 
is present. 

The route to the usual GPE is to assume that the full 
field operator can be replaced by a wave function ■(/'(x) — 
i.e. that all quantum fluctuations can be neglected. We 
proceed instead by defining a projection operator V such 
that 



P'i'(x) = «k</'k(x). 



(3) 



where the region C is determined by the requirement 
that (aj^ak) ^ 1, and the set {(/>k} defines some basis 
in which the field operator is approximately diagonal at 
the boundary of C. For these modes, the quantum fluctu- 
ation part of the projected field operator can be ignored, 
and so we replace at Ck and write 



Y Ck(/>k(x). 



(4) 



Defining the operator Q — 1 ~ V and Q^(x) = ?)(x), 
operating on Eq. (|^) with V and taking the mean value 
results in what we call the finite temperature GPE 



ih 



dt 



7?,pV(x)+[/o^{|^/'(x)|2^(x)} 

C/o75{2|^(x)p(,)(x))+^(x)2(r)t(x))} 

C/oP {r (x) (ry(x)r7(x)> + 2^(x) (t?^ (x)7}(x)) } 

C/oP{(f/t(x)7}(x)r7(x)>}, (5) 



This equation describes the full dynamics of the coherent 
region and its coupling to an effective heat bath described 
by ?y(x). In general, the non-equilibrium evolution de- 
pends on the coupling between these two regions and the 
exchange of energy and particles that this allows. The 
FTGPE must be complemented by an equation of mo- 
tion for 77(x) and in principle this can be obtained using 
a form of quantum kinetic theory. 

The only approximation that has been made in the 
derivation of the FTGPE is that the modes represented 



by V'(x) must satisfy the criterion of classicality, that is 
Nk ^ 1. The FTGPE is a non-perturbative equation, 
and therefore we expect that it will be valid in the re- 
gion of the phase transition as long as only the highly 
occupied modes are treated. There is perhaps a misper- 
ception in the BEG community that the GPE is only 
valid at T = 0. However, it is well known that close to 
the phase transition a classical description of the longer 
length scales involved is completely appropriate. This is 
exactly what the GPE describes, and in fact it has been 
used as a model of phase transitions in other areas of con- 
densed matter physics. Indeed our model has the same 
energy functional for these modes as used in the classi- 
cal renormalization group theory of the superfluid phase 
transition. It therefore seems reasonable to expect that 
the same approximations are valid in this case. 

The physical processes described by the various terms 
of Eq. (|) are discussed in detail in Ref. . In this paper 
we concentrate on a simplification of the FTGPE which 
is effectively a model of a restricted system. This allows 
us to demonstrate some of the properties of the GPE 
without having to solve the more complicated equation. 



III. THE PROJECTED GPE 



In this paper, we wish to show that the GPE alone 
can describe evolution of general configurations of the 
coherent region C towards an equilibrium that can be 
parameterised by a temperature. We therefore ignore all 
terms involving ■q{yi) in Eq. (||) and concentrate on the 
first line 



iTi 



aV^(x) 
dt 



7?,pV^(x)-HC/o7^{|V(x)PV'(x)}, (6) 



which we call the projected Gross-Pitaevskii equation 
(PGPE). Although Eq. (|6|) is completely reversible, it 
is well known that deterministic nonlinear systems with 
only a few degrees of freedom exhibit chaotic, and hence 
ergodic behaviour p3| . If many modes are occupied, the 
PGPE contains many degrees of freedom and it is there- 
fore reasonable to expect it to evolve to equilibrium (ex- 
cept for specially chosen initial conditions such as eigen- 
state solutions). 

The projected GPE describes a microcanonical system. 
However, if the region C is large, then its fluctuations in 
energy and particle number in the grand canonical en- 
semble would be small. Hence we expect the final equi- 
librium state of the projected GPE to be similar to that of 
the finite temperature GPE coupled to a bath 77 (x) with 
the appropriate chemical potential and temperature. The 
detailed non-equilibrium dynamics of the system will de- 
pend on the exchange of energy and particles between C 
and the bath — however, we leave the coupling of 
and J7(x) to be addressed in future work. 
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A. The projector 

The spatial representation of the projection operation 
is written 



7'{F(x)} 



^fc(x) / d^x' 0Ux')i^(x'), (7) 



and this operation must be carried out numerically every 
time we calculate the nonlinear term in the PGPE. This 
is a very time consuming operation in general, taking 
many times longer than calculating itself. 

The operation is much simpler numerically if we use a 
plane-wave basis in our projector 



exp(ik • x) 



(8) 



where V is the volume of our system. In this case 
Eq. (|^) becomes simply the application of a forward 
Fourier transformation to our function i^(x), followed by 
an inverse Fourier transformation that includes only the 
modes in the coherent region. Thus our numerical pro- 
cedure is 



'P{F(x)} = IFFT<^ P(k) X FFT [F(x)] 



(9) 



where FFT and IFFT refer to the forward and inverse 
fast Fourier transform operations respectively, and P(k) 
is the representation of the projector V in Fourier space. 
There are very efficient routines available to carry out 
FFTs, and so we find that it is extremely advantageous 
numerically to define our projector in the plane- wave ba- 
sis. 



changes with the condensate fraction. In general the con- 
densate fraction must be determined by diagonalisation, 
which can be a very time consuming procedure pO[. 



IV. SIMULATIONS 

We have performed simulations for a fully three- 
dimensional homogeneous Bose gas with periodic bound- 
ary conditions. The dimensionless equation we compute 
is 

= -VV^i) + CniT'il^WlVW}, (10) 

where the normalisation of the wave function has been 
defined to be 



J d3x|V'(x)P 1. 



The nonlinear constant is 



_2mNUo 



(11) 



(12) 



where N is the total number of particles in the system, 
and L is the period. Our dimensionless parameters are 
X — x/L, wave vector k = kL, energy e = s/el, and 
time r = SLt/h, with sl = / {2mL'^) . 



A. Parameters 

The two parameters that determine all properties of 
the system are the projector P and the nonlinear con- 
stant Cnl. 



B. Implications 

For any non-periodic trapping potential, the use of a 
plane- wave basis is at odds with our requirement that the 
basis must approximately diagonalise ipijn) at the bound- 
ary of the region C. In fact, it may not even satisfy this 
requirement for a periodic potential if the boundary of 
the coherent region occurs at a low enough energy. 

If we consider a homogeneous system, however, the 
plane- wave basis will always satisfy our requirements. In 
this case the effect of a condensate on the excitations of 
the system is simply to mix modes of momenta p and —p. 
Thus even if V'(x) is not diagonalised at the boundary of 
C, we can still apply the projector cleanly in Fourier 
space. For these reasons, the simulations that we present 
in this paper are for the homogeneous Bose gas. We 
intend to address the issue of projectors for the trapped 
Bose gas in future work. 

A direct advantage of simulating the homogeneous sys- 
tem is that the condensate occupation is readily identified 
as the k = component of the wave function. This is in 
contrast to the trapped case, where the condensate mode 



1. Projector V 

We have chosen a projection operator such that all 
modes included in the simulations have |k| < 15 x 2tt/ L, 
which enables us to use a computationally efficient nu- 
merical grid of 32 X 32 X 32 points. This means that 13997 
modes are included in the system. 

Grid size and aliasing The nonlinear term of the GPE 
can generate momentum components up to three times 
larger than those which exist in the original wave func- 
tion. Thus it would seem that calculating the term 
|?/'(x)p?/'(x) on a grid only slightly larger than the pro- 
jector would cause problems with aliasing. The correct 
procedure would be instead to calculate this term on a 
grid size of 96 x 96 x 96 points before performing the 
projection operation. 

To check the effect of grid size we have performed simu- 
lations where the nonlinear term was calculated on grids 
of size 32, 64, and 96 points, and found that there is 
no difference in the equilibrium properties of the system. 
The detailed dependence of the condensate population 
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during evolution is different in detail for each size grid, 
but follows the same average curve. The same behaviour 
is observed when adjusting the accuracy parameter of our 
adaptive step size algorithm for evolving the GPE. 

We attribute this behaviour to the deterministic chaos 
exhibited by the system. Any small numerical error is 
eventually magnified such that the system follows a quite 
different microscopic path through phase space, although 
the resulting macroscopic (average) properties are unaf- 
fected. 



2. Nonlinearity Cni 

We note that the choice of the nonlinear constant de- 
termines only the ratio of NUq/L. This means that for 
a given value of Cni, we can choose the parameters N, 
Uq and L such that our condition A'k = A^jckp ^ 1 is 
always satisfied for a given physical situation. 

We have performed three series of simulations with 
nonlinearities of Cni = 500, 2000, and 10000. The high- 
est value of Cni was chosen such that all the states con- 
tained in the calculation are phonon-like for a large con- 
densate fraction. The boundary between phonon-like and 
particle-like states for the homogeneous gas is 



2m 



= n-oUo, 



(13) 



where we have defined A^o to be the condensate number 
within the volume L^, and thus no = Nq/L^ is the con- 
densate density. Converting Eq. (|l^) to dimensionless 
units we find that 



ko = -\/Cnl — , 



(14) 



and therefore for a condensate fraction oi Nq/N = 1 we 
have 

C„i = 10000 ^ fco » 15.9 X 2n, 
C„i = 2000 fco « 7.12 X 2n, 
C„i = 500 fco « 3.56 X 2n. 

We find that computations with smaller values of Cni 
take comparatively longer to reach equilibrium. This is 
because the equilibration rate is approximately propor- 
tional to C^j, whereas the minimum time-step required 
for a given accuracy in the numerical integration of the 
PGPE only increases slowly with decreasing Cni. 

To give an indication of how these dimensionless pa- 
rameters compare to experimental setups, for Cni = 
10000 we can choose ^"^Rb atoms with A = 1.8 x 10^ 
and L « 26 /zm to give a number density of about 10^"* 
cm""^ — similar to current experiments on BEC in traps. 



B. Initial wave functions 

We begin our simulations with strongly non- 
equilibrium wave functions with a chosen total energy 



E. We construct these by populating the amplitudes of 
the wave function components Ck in the expansion 



V'(x,0) 



(15) 



The populations |ckp are chosen such that the distribu- 
tion is as flat as possible, while the phases of the ampli- 
tudes are chosen at random 

The total energy i5 is a constraint on the distribu- 
tion of amplitudes. The energy of a pure condensate is 
Eq = Cni/2, all of this being due to interactions — the 
kinetic energy is zero. To have a wave function with an 
energy not much larger than Cni/2, the occupations of 
the k = state and the fc = 27r states cannot be equal. 
(We use the notation k = |k|.) Therefore, for the lowest 
energy simulations the initial condensate population is 
necessarily larger than the excited state populations. 

To ensure that the initial wave functions are sufficiently 
randomised, we enforce the condition that all 123 states 
with fc < 3 x 27r must have some initial population, while 
all other components may be unoccupied. For low en- 
ergies, when this distribution including the condensate 
cannot be totally flat, we keep the populations of the 
components with 1 < fc/27r < 3 equal, and adjust the 
condensate population such that the wave function has 
the energy we require. An example of this situation is 
shown in Fig. |l|(a) for the E = 7000 initial wave function 
in the Cni = 10000 simulation series. 

For simulations with a sufficiently high total energy E 
that the inner 123 components may have equal popula- 
tion, we continue to add further shells of higher k to our 
wave function. The amplitudes of the inner components 
are readjusted to maintain the required normalisation. 
This causes the energy of the system to increase mono- 
tonically with each new shell until we find two wave func- 
tions that bound the energy we are looking for, differing 
only in their outermost shell. We then adjust the pop- 
ulation of the outermost shell downwards until we reach 
the required energy. 

This procedure is necessary due to the nonlinearity of 
the problem. In the case of the ideal gas (Cni — 0), we can 
calculate the kinetic energy (and hence the total energy) 
of the wave function simply by knowing the distribution 
of |ckp, via 

+-2 p 

Eun = -— / d3x^*(x)VV(x), 



2m 



(16) 



However, for Cni > we must also add the interaction 
energy of the wave function to the total energy. This is 



— El \ " * * X 

2 / y CpCqCjnCnOp+q— m— nj 
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(17) 
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FIG. 1: Two dimensional slices of wave functions through the 
kz — plane in momentum space for the Cni = 10000, E — 
7000 simulations, (a) Base 10 logarithm of the fc-space wave 
function at r = (b) Base 10 logarithm of the fc-space wave 
function at r = 0.2 once the system has reached equilibrium. 



and depends nontrivially on the {ck}. 

Further images of initial and final state wave functions 
are shown in Fig. |l| in A;-space, and Fig. || in real space. 



C. Evolution 

The PGPE is evolved in the interaction picture, us- 
ing a fourth-order Runge-Kutta method with adaptive 
step size determined by estimating the fifth-order trun- 
cation error. The acceptable relative truncation error 
was set to be 10^^*^ for all components with an occupa- 
tion of > 10"'*A*'o/A^. This resulted in typical time steps 
as presented in Table |, which could be integrated in a 
reasonable time on a modern workstation. 

We evolve the initial wave functions for at least twice as 
long as it takes for the system to reach equilibrium, based 
on the observation of the behaviour of the condensate 
fraction (see Sec. The time period for each value of 
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FIG. 2: Two dimensional slices of wave functions near the 
z = plane in real space for the Cni = 10000, E = 7000 simu- 
lation, (a) Base 10 logarithm of the real-space wave function 
at r = (b) Base 10 logarithm of the real-space wave function 
at r = 0.2 once the system has reached equilibrium. 



Cni 


Min. time step 


Max. time step 


Length of 


(io-«) 


(10-6) 


evolution r 


500 


4 


6 


2.0 


2000 


1.6 


4.4 


0.4 


10000 


0.45 


1.2 


0.2 



TABLE I: The typical minimum and maximum time steps for 
the simulations. The minimum is for high energy simulations, 
and the maximum for low energy. 



Cni is also given in Tabic |. Thus the longest of these 
simulations required ~ 5 x 10^ time steps. 



V. EVIDENCE FOR EQUILIBRIUM 

Although the PGPE is completely reversible, the final 
state wave functions displayed in Figs. ^ and ^ indicate 
that the simulations have evolved the system to an ap- 
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parent equilibrium state. The fc-space distributions have 
evolved from initially being flat to a form that is peaked 
at the centre, and tails away towards the edges. Also, 
there is a smoothing out of both the phase and density 
profiles of the real-space wave function. After a certain 
time of evolution Toq, the plots for the wave functions 
appear to be isomorphic for r > Teq. 

We would like to note that the equilibrium proper- 
ties depend only on the total energy and momentum of 
the initial wave function — they are independent of the 
shape of the initial distribution in fc-space. We have per- 
formed simulations with non-spherical initial wave func- 
tions, and found that they evolve to a spherical equilib- 
rium state. Also, as the GPE conserves momentum, for 
the condensate to form in the k = mode the initial 
distribution must have zero total momentum. We have 
performed simulations where the initial distribution had 
a finite momentum, and observed the condensate to form 
in a non-zero momentum state as expected. 

To determine the properties of the system at equilib- 
rium, in theory we should carry out many different sim- 
ulations each with the same initial populations but with 
different choices of the initial phases, and then take the 
ensemble average. However, this is an extremely large 
computational task. Instead, we assume the ergodic the- 
orem applies, such that the time average over the evo- 
lution of a single system at equilibrium is equivalent to 
the ensemble average over many different systems. We 
therefore perform a time-average over the last 50 wave 
functions saved, all with r > Teq. 



1. Condensate occupation 

Strong evidence that the simulations have reached 
equilibrium is provided by the time dependence of the 
condensate population. For all simulations this settles 
down to an average value (dependent on the energy E) 
that fluctuates by a small amount. The initial time evo- 
lution of the condensate fraction for five different energies 
with C„i = 10000 is shown in Fig. |. 

The average condensate occupation in equilibrium for 
all simulations for the Cni — 10000 case are presented in 
Fig. ^(a). The fluctuations of the condensate population 
are indicated by the (barely visible) vertical lines at each 
point, and these are largest for the E = 9000 simulation. 
For comparison, the corresponding curve for the ideal 
gas is plotted in Fig. ^(b). We can see that for Cni — 
the curve is linear up to the transition point, but the 
Cni = 10000 curve displays a distinct bulge in this region. 
The shape of the corresponding curves for Cni = 500 and 
2000 faU in between the C„i = and 10000 cases. 



2. Particle distribution 

Further evidence of equilibrium is provided by the dis- 
tribution of the particles in momentum space. Rather 
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FIG. 3: Plot of the initial time evolution of No{t)/N for four 
different simulation energies with Cni = 10000. From top to 
bottom: E = 5500, 7000, 8500, 9250, 10000. The simulations 
were run until r = 0.2. Other values of the nonlinearity give 
qualitatively similar results. 



than using the plane- wave basis, we transform the wave 
functions into the quasiparticle basis of quadratic Bogoli- 
ubov theory. For the homogeneous gas, this theory can 
be solved analytically and we can write the quasiparticle 
amplitude 6k as 



where 



-ak 



,.2 ' 



and ak is given by 



vl- 



(18) 
(19) 

(20) 



In this last equation, the dimensionless wave vector 
is given hy yk = k/ko with fco as defined in Eq. (|lij). 
The normalisation condition — = 1 is automatically 
satisfied by Eq. ([l9|). From Eq. (14) we can see that the 
sole parameters of the transformation are the condensate 
fraction {No)/N, and the nonlinear constant Cni. 

We time average the populations of the quasiparti- 
cles states Ny^/N = jfokP as was described in Sec. ^ to 
give (A'k) /N, and finally average over angle so that we 
can produce a one dimensional plot of {Nk)/N. This 
distribution for four different simulation energies and 
C„i = 10000 is shown in Fig. |. 

We can see that the shape of the curves is surprisingly 
smooth for each energy, suggesting that the system is in 
equilibrium. The plot of the distribution for any individ- 
ual wave function is scattered about the average. 

We have also determined the fluctuations of the pop- 
ulation of the quasiparticle modes. The grand canonical 
ensemble for the Bose gas predicts the relationship 



(21) 
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FIG. 4: (a) Condensate fraction plotted against total energy 
after each individual simulation has reached equilibrium for 
Cni = 10000. The barely discernible vertical lines on each 
point indicate the magnitude of the fluctuations, (b) The 
curve for the same system, but calculated for the ideal gas. 



for fc 7^ 0, which in the classical limit (Nk) ^ 1 gives 

(ANk) ~ (Nk), (22) 

This is indeed the behaviour that we observe. Although 
we are evolving a microcanonical system, in this case 
there are such a large number of modes that the remain- 
der of the system acts as a bath for any individual mode 
and the result still applies. 



VI. QUANTITATIVE ANALYSIS OF THE 
DISTRIBUTIONS 

While the data presented in Sec. ^ indicates that the 
PGPE is evolving the system to equilibrium, as yet we 
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FIG. 5: Plots of the equilibrium Bogoliubov quasiparticle dis- 
tributions averaged over time and angle for four different to- 
tal energies. Squares E — 6000, crosses E = 7500, circles 
E — 9000, dots E — 11000. The mean condensate occupation 
for the first three distributions is off axis. 



have presented no quantitative evidence. To demonstrate 
conclusively that equilibrium has been reached, we need 
to be able to assign a temperature to the simulations. In 
this section we measure a temperature for a given simu- 
lation by comparing the distribution function of the nu- 
merical simulations against a predicted energy spectrum. 



A. Expected equilibrium distribution 

The GPE is the high occupation (classical) limit of 
the full equation for the Bose field operator, Eq. (|^). 
Therefore, in equilibrium we expect the mean occupation 
of mode k to be the classical limit of the Bose-Einstein 
distribution — i.e. the equipartition relation 



£k- 



(23) 



where k labels the eigenstates of the system. In general 
these will be some type of quasiparticle mode. Manipu- 
lating Eq. (|2^) , we find that 



£k 



keT 



+ ti. 



(24) 



The equilibrium condensate occupation according to the 
equipartition relation will be given by Eq. (p3) with 
(Nk) (Nq) and A (the condensate eigenvalue). 

From this expression we can solve for the chemical po- 
tential 



A 



keT 
(No)' 



(25) 
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Substituting this result into Eq. 
dimcnsionless units wc find 

4- A ( N 



and converting to 



N 



\{Nk) (No) 



(26) 



where T ~ ksT/ (Nel) is the dimensionless temperature. 

Once equihbrium has been reached for a single simu- 
lation, we make use of Eq. ( p6| ) to measure the quantity 
T. Decomposing the wave functions in some basis and 
time-averaging the populations determines (Nk) IN as a 
function of the variable k, as is plotted in Fig. This 
completely specifies the RHS of Eq. ( |26| ) and it remains 
to determine the quantities on the LHS. 

In this section we consider three different methods of 
either predicting or measuring the function ik — A. If 
the basis we have used for our decomposition is a good 
one, and our prediction for efe — A is correct, then this 
curve will have the same shape as the RHS of Eq. (26). 



The constant of proportionally determined by a fitting 
procedure will then give the temperature T. 

Before we describe our methods and results, we would 
like to note that the quantity we refer to throughout the 
remainder of this paper as the temperature is the vari- 
able T as determined by the numerical fitting procedures 
described above. We have not yet established that this 
is the true temperature as defined by thermal equilib- 
rium with a heat reservoir. However, we believe that if 
we were to solve the FTGPE with 77(x) acting as a heat 
bath, then the temperature determined in the coherent 
region via this method would agree with the bath tem- 
perature. 



B. Method 1 : Bogoliubov theory 

In the limit of large condensate fraction {No)/N ^ 1, 
we expect the Bogoliubov transformation to provide a 
good basis. For the homogeneous case the dispersion 
relation is analytic, and is given by 



2m 



-,1/2 



{chkf 



(27) 
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FIG. 6: Fits of the simulation quasiparticle population data 
to the Bogoliubov dispersion relation for two cases. For both 
graphs the solid line is the Bogoliubov curve, while the dashed 
line is the ideal gas dispersion relation. The temperature is 
determined by a least-squares fit to the plot of {N/{Nk) — 
N/{No)), which is shown as the dots, (a) Cni = 500, E = 
500 and {No)/N = 0.929, with a best fit temperature from 
Bogoliubov theory of f = 0.0175. (b) Cni = 10000, E = 
5250 and {No)/N = 0.957, with a best fit temperature from 
Bogoliubov theory of f = 0.018. 



where c = (ngiJo/m)-^/^ is the speed of sound and Sk is 
the absolute energy of a mode with wave vector k. In our 
dimensionless units this becomes 

ek-~X=(k^ + 2cJ-^py\ (28) 

The condensate fraction is determined from the numer- 
ical results, and so when the Bogoliubov dispersion re- 
lation is valid we can determine a temperature for the 
simulations by substituting Eq. ( p8| ) in Eq. ([2^). 



Results 

We have carried out this analysis for all the simulation 
data. For the Cni = 500 case, the measured distributions 
are in excellent agreement with the Bogoliubov dispersion 
relation for all energies, and we have been able to extract 
the corresponding temperature for each simulation. 

However, this is not the case for the more strongly 
interacting systems. For Cni = 2000, the Bogoliubov re- 
lation is a good fit only for simulations with E < 2000 
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{{No)/N > 0.75), or for energies above the BEC tran- 
sition point. For the C^i — 10000 case, good agree- 
ment is found only for the lowest energy simulation with 
E = 5250 and {No)/N « 0.96. Sample fits of the sim- 
ulation data to the Bogoliubov dispersion relation are 
shown in Fig. ^ for cases where the agreement is good. 
(An example of this procedure for where the Bogoliubov 
spectrum is not appropriate is given in Fig. P). 

The reason for the limited range of agreement is be- 
cause the Bogoliubov transformation diagonalises only a 
quadratic approximation to the full Hamiltonian. It ne- 
glects terms that are cubic and quartic in non-condensate 
operators, assuming that they are small (these are dis- 
cussed in detail below). This is a good approximation for 
the Cni — 500 simulations — at large condensate fraction 
the dispersion relation is only slightly shifted from the 
non-interacting relation ik = , and at smaller conden- 
sate fractions the difference is negligible. Hence we can 
fit a temperature up to and above the BEC transition. 

For the Cni — 2000 case the higher order terms be- 
come important above E = 2000, and for the strongest 
interaction strength of Cni — 10000, they are important 
for all but the lowest energy simulation we consider. For 
the higher energy simulations the shape of Eq. ( p6| ) no 
longer agrees with Eq. (|2^), and we must use a more 
sophisticated theory to predict the dispersion relation. 

Above the transition point, however, there is no con- 
densate and the ideal gas dispersion relation is a reason- 
able description of the system. 



C. Method 2 : Second order theory 

As the occupation of the quasiparticle modes becomes 
significant at large interaction strengths, the cubic and 
quartic terms of the many-body Hamiltonian that are 



neglected in the Bogoliubov transformation become im- 
portant. In Ref. 1^ Morgan develops a consistent exten- 
sion of the Bogoliubov theory to second order that leads 
to a gapless excitation spectrum. This theory treats the 
cubic and quartic terms of the Hamiltonian using pertur- 
bation theory in the Bogoliubov quasiparticle basis. This 
results in energy-shifts of the excitations away from the 
Bogoliubov predictions of Eq. (27). 

Expressions for the energy-shifts of the excitations are 
given in Sec. 6.2 of Ref. They have the form 



Aik = AEsik) + AEi{k) + AExik), 



(29) 



where AE^ik) [Ai?4(fc)] is the shift in energy of a quasi- 
particle in mode k due to the cubic [quartic] Hamiltonian, 
and AE\{k) describes the shift due to the change in the 
condensate eigenvalue. In the high-occupation limit we 
find 

AEiik) + AE^ik) = -CniAi^li^, (30) 



1 



at 



where k is the dimensionlcss anomalous average, defined 
by 



E 



{Nk + N-k)ak 



N{1 



(31) 



The expression for Ai?3(fc) is derived from second- 
order perturbation theory, and is rather complicated. We 
have 



A^3(fc) 



where 



— 2Cnl 

1-a?. 



[A£;^(fc) -I- AEl(k) -f AEl{k)], (32) 



AE'Hk^ = [N', + Nj){l -aj- + a^ak + ajak ~ aiajakY /oo^ 
^ Y N{z, + z,-Zk){l-aiY{l-a,Y ' ^ ^ 

' y N{z, + Zj+Zk){l-a.,f{l-ajY ' ^ ' 

AE'^lk) = V ~ ^J^^-^ ~ ^ "fc + + Q^'^fc - "^"jQfc)^ ^ 

^ Y N{z,^ Zj+Zk){l-a,)^l-a,y ' ^ ^ 

in which i = k — j , and 

Zk = yk(2 + ylf/^ ^ ek{c.^-^y\ 

is another form of the dimensionless energy of mode k, with tjk ~ k/ko as earlier. 

I 



(36) 



1. Calculation of energy shifts the shifts for our simulations. The first procedure is to 

The numerical calculation of the energy shifts is not 
a trivial task, and have used two methods to determine 
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FIG. 7: The shifts to the Bogoliubov quasiparticle energies for 
two different simulations. The soUd thin curves are calculated 
via the first method described in the text using population 
data extracted from the simulations, and are hence somewhat 
noisy. The thick grey curves use the second method, assum- 
ing equilibrium populations given by Bogoliubov theory and 
calculated by numerical integration. The lower curves are 
for the Cni = 2000, E = 4000 simulation, and appear to be 
approximately gap less as 0. The upper curves are for 

the Cni = 10000, E = 6000 simulation, and exhibit a gap as 
fc ^ 0. 



calculate the shifts directly using the population data 
from the simulations. We therefore 

1. Calculate the quasiparticle populations for the 
last 50 wave functions of our simulation based on 
a condensate population (iVo), and then average 
these over time. 

2. Calculate the energy shifts for mode k using these 
populations as the input. 

3. Average the shifts over angle to give a one- 
dimensional function of k. 

This results in plots of the energy shifts that are some- 
what scattered due to the finite size of the system. The 
expressions for the shifts Eqs. (|32-35) contain poles when 
energy matches occur, and hence the numerical calcula- 
tion is performed using an imaginary part in the denomi- 
nator. The size of this imaginary part does not affect the 
shape of the curve in the limit that it is small, but it does 
affect the amount of scatter in the shifts. We have per- 
formed sample calculations allowing L to increase while 
keeping other parameters of the system constant, and 
this makes the curve smoother. 

The second procedure only makes use of the con- 
densate fraction and the total number of quasiparticles, 
rather than the population of the individual levels. By 
assuming the Bogoliubov spectrum is a good estimate 
of the energies (which must be true for the perturbation 



theory to be valid) , we can estimate the temperature Tost 
using the normalisation constraint on the populations 



E 



N 



N 



fe>0 



£k - A 



(37) 



where we have used the approximation fl — X that is 
valid when there is a condensate present. The LHS as 
well as the value of {No)/N are determined by the simu- 
lations, and the Bogoliubov relation Eq. ( p8| ) is used for 
the energies. 

Once the estimated temperature Tost is determined, we 
use the equipartition Bogoliubov relation for the popula- 
tions in Eqs. (^3|-^5|), and then approximate the sums by 
numerical integration to calculate the shifts to the levels. 
We find that this gives curves that agree on average with 
those calculated using the first method, but are much 
smoother. A comparison of the two methods is given in 
Fig. 0. 



2. Results 

For the Cni = 2000 simulations, the quasiparticle pop- 
ulations extracted from the simulations are in much bet- 
ter agreement with the energy spectrums from the second 
order theory than with those from ordinary Bogoliubov 
theory. We find that most of the measured distributions 
for the Cni = 2000 case are well described by the second 
order theory. Sample results are presented in Fig. ^(a). 

However, this is not the case for the Cni = 10000 sim- 
ulations. In fact we find that the energy spectrum is 
shifted in the opposite direction to that inferred from the 
simulations, and that there is a energy gap for fc 0. 
The reasons for this are discussed below. 



3. Breakdown of perturbation ttieory 

The validity of the second order theory is constrained 
by the requirement H 



(38) 



where no is the condensate density. This corresponds in 
our dimcnsionless units to 



T 



Cni 



(87r)3/2 \{Nq)/N 



1/2 



< 1. 



(39) 



For the results of Fig. | with Cni = 2000, E = 4000 
this parameter is 0.14 and so we are beginning to probe 
the boundary of validity of the theory. At higher E the 
shifts become of the order of the unperturbed energies, 
and hence the results are unreliable. In this region even 
higher order terms are important, and the second order 
theory can no longer be expected to give good results. 
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FIG. 8: Fits of the simulation quasiparticle population data 
to dispersion relations. The dots are a plot of {N/{Nk} ~ 
N/{No)), the solid curve is for dispersion relation predicted 
by second order theory, and the dashed curve is the dispersion 
relation predicted by Bogoliubov theory, (a) Cni = 2000, 
E = 4000, and {No)/N = 0.279. Second order theory gives a 
good fit to the numerical results with a best fit temperature of 
f = 0.201. (b) Cni = 10000, E = 6000, and {No)/N = 0.841. 
The shape of the second order theory dispersion relation does 
not agree with the population data from the simulation, and 
the gap is apparent as fc ^ 0. The grey curve plot s the 
energies as determined by the inethod described in Sec. VI D 
with a best fit temperature of T = 0.0726. 



From our calculations it seems that this parameter should 
be < 0.2 for the theory to be valid. 

We would like to emphasize, however, that the GPE 
suffers no such limitations. It is non-perturbative and 
thus we expect that it will be valid all the way through 
the transition region as long as the high occupation num- 
ber condition is satisfied. 



4- Gaplessness in a finite system 

In the course of this work it has become apparent that 
while the second order theory is gapless for infinite sys- 
tems, this is not the case for systems such as ours with 
a finite momentum cutoff. The individual terms in the 
perturbation expansion given by Eqs. ( ^0| ) and con- 
tain contributions that are proportional to 1/fc (infrared 
divergent) and a constant (gap) in the low k limit. For 
a homogeneous system these terms cancel exactly when 
the upper limit of the integrals is infinite, and this leaves 
a gapless spectrum ^, 0. However, in a system with a 
momentum cutoff these terms do not exactly cancel, with 
the result that there is a gap in the predicted excitation 
spectrum as fc — > 0. 

Briefly, this gap arises because the energy shifts 
AEi{k) + AExik) of Eq. (|^) only involve the quan- 
tity K. This is obtained from Eq. ( |3l] ) via a sum over 
all states below the cutoff where the summand depends 
only on a single wavevector. In contrast the shift Aii^3(fc) 
of Eq. (^2|) involves a sum over states where the sum- 
mand depends on two wavevectors i and j (related by 
momentum conservation) both of which must be below 
the cutoff. This difference in the restrictions on the sum- 
mations leads to a lack of complete cancellation in the 
corresponding shifts at low energy and the appearance of 
a gap in the excitation spectrum. 

For the homogeneous gas, we can calculate the size of 
the gap predicted by the second order theory analytically. 
Replacing the summations by integrations, we find that 
the leading order contribution to the energy shift in the 
limit fc — !■ is 



Aefc 



noUo 



(rioa ) 



3U/2 



1/2 



. (40) 

where y — k/kQ as before and yc — fcc/^o where k^ is the 
momentum cutoff. In the limit k ^ we have cx yk 
and so Ae^ tends to a constant (the gap). The size of 
the gap tends to zero as the momentum cutoff j/c tends to 
infinity but otherwise it is finite. We stress that Eq. ( |40| ) 
is only the low k limit of the exact result. For our simu- 
lations there is a minimum wavevector in the problem so 
it is possible for the terms of order yk to be larger than 
the gap contribution given above. This is the case for the 
simulations with Cni ^ 2000 where j/c is reasonably large 
and the gap is therefore small. 

The result of Eq. (|4^) contains the small parameter 
that controls the validity of the second order theory in 
the usual case where there is no momentum cutoff [c.f. 
Eq. (^)]. However, the result also depends exphcitly on 
the cutoff kc so in this case there is a second parameter in 
the theory. For perturbation theory to be valid we require 
that the predicted energy shifts are small compared to 
the unperturbed energies, i.e. that Ae^/efc <C 1. We 
therefore obtain a second criterion for the validity of the 
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second order theory which is 



(Ml 

\noUo 



{noa 



3N1/2 



1/2 



1 



71"/ Vki'i + yl) 



< 1. 



(41) 



This resuh should hold for all momenta in the simula- 
tions, and in particular for the smallest value of yk- For 
the Cni = 2000, E = 4000 simulations the left hand 
side is 0.04 for k = 2tt. In this case the gap is negligi- 
ble and the dominant contribution to the energy shifts 
comes from the terms of order yk in Eq. (^ ) . The small 



parameter of the theory is therefore given by Eq. (38). 
However, for the C^i = 10000, E = 5250 simulations the 
left hand side is of order 0.12. In this case the gap is 
not negligible and we cannot use second order theory to 
define a temperature. 

This result is somewhat surprising since, even for 
a condensate fraction of 80% the small parameter of 
Eq. ( |3^ ) is of order 0.07, and it does seem reasonable 
to expect that perturbation theory should be applicable. 
This does not appear to be the case, however, and we 
have so far been unable to determine the root cause of 
this problem. It is worth noting that the numerical simu- 
lations themselves have no difhculties in this regime and 
do not predict a gap at low momentum. This is because 
the GPE is non-perturbative and indeed this is one of 
the main reasons for using it to study the properties of 
Bose condensed systems at finite temperature. 

The disagreement between the second order theory 
and the numerical simulations is illustrated in Fig. ^(b), 
where it can be seen that even despite the gap, the shifts 
the theory predicts are in the wrong direction in compar- 
ision with the simulations. 



D. Method 3 : Non-perturbative determination of 
the temperature 

The failure of second order theory for the Cni = 10000 
simulations caused us to investigate other possible meth- 
ods of determining the temperature once the system was 
in equilibrium. This has lead to what seems to be a 
method of determining the temperature that does not 
rely on perturbation theory, and we describe it here. 

We found earlier that the Bogoliubov spectrum gave 
a good prediction of the populations of the quasiparticle 
levels for the lowest energy simulation in the Cni — 10000 
series with E = 5250. Therefore it seems reasonable 
that the Bogoliubov basis should remain a good one for 
perturbation theory for the next simulation with E = 
5500, even though the second order theory cannot be 
used to calculate the energy shifts. 

Therefore we attempted another method to determine 
the absolute energy of each quasiparticle level. If we 
are using a good basis, then on a short time scale the 
quasiparticles should be independent, with amplitudes 
evolving according to 



&k(0 = &k(io) exp(-i£k(i - tQ)/h). 



(42) 



Thus by measuring the gradient of the phase of each 
quasiparticle we can determine its energy. 

To determine the energy spectrum for a single simula- 
tion, our numerical procedure was as follows. 

1. Take the last 50 wave functions saved for a simu- 
lation once it has reached equilibrium, and evolve 
each of these individually for a very short period. 
One hundred wave functions are saved for each of 
the 50 simulations. 

2. Transform the wave functions into the quasiparticle 
basis, and measure the energy of each quasiparti- 
cle, determined by a linear fit to the phase of each 
amplitude for all 50 simulations. 

3. Average over all 50 energy spectrums to give a sin- 
gle three dimensional spectrum. 

4. Finally, average over angle to give a one dimen- 
sional energy spectrum. 

This gives us a dispersion relation — A which can 
then be compared to a plot of {N/(Nk) - N/(Nq)). If 
the shapes of the curves agree, then a temperature can 
be determined via Eq. ( p6| ) as in the earlier sections. 

We first tested this procedure on the Cni = 2000 sim- 
ulation series, and found that this method was in good 
agreement with the second order theory calculations, the 
two approaches assigning the same temperature to the 
various simulations. 

We then moved onto the Cii = 10000 simulations. We 
found the surprising result that not only did the shape of 
the plots of the curves for efc/fcsTfit and (1/A^fc — I/Nq) 
agree for the lower energy simulations where the parame- 
ter of Eq. (|3|) was small, it also agreed when it was of the 
order of, and greater than one. This was unexpected, as 
it would seem likely that near the phase transition when 
interactions are strong that the Bogoliubov quasiparticle 
basis would no longer be sufhciently good for this method 
to be accurate. An example of the energy spectrum and 
its fit to the population data is shown in Fig. ^(b). 

As a further test we carried out the same procedure de- 
scribed above, but using the plane wave basis rather than 
the Bogoliubov quasiparticle basis. Intuitively it would 
seem that this would no longer work — but we found that 
not only did it give the same temperatures as the quasi- 
particle basis for the Cni = 10000 simulations, it also 
agreed with the temperatures determined using second 
order theory for the Cni = 2000 simulations. 



VII. TEMPERATURE DEPENDENCE 

Using the three methods described in the previous sec- 
tion, we have been able to measure an equilibrium tem- 
perature for all simulations in this paper. We are con- 
fident of the results determined from both Bogoliubov 
and second order theory; however unfortunately we do 
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FIG. 9: Condensate fraction versus temperature for the 
PGPE system with fc < 15 x 2tt/L for four different interac- 
tion strengths. The open circles are for Cni = 10000, crosses 
for Cni = 2000, sohd dots for Cni = 500, and the soUd hne is 
for the ideal gas. The shift in the transition temperature is 
positive with increasing interaction strength Cni. 



gas dispersion relation, and therefore the transition tem- 
perature will not be greatly shifted over a wide range of 
nonlinearities. 

There has been some discussion recently in the liter- 
ature about the shift in the transition temperature for 
the homogeneous interacting Bose gas, with some au- 
thors even disagreeing in the direction of the shift (e.g. 
see Ref. and references within). For the PGPE sys- 
tem described in this paper by Eq. (^), we can see from 
Fig. H that the shift is positive, although as yet we have 
made no effort to quantify this. This would require many 
more simulations to be run, especially in the transition 
region, and for the temperatures to be determined more 
accurately. 

It seems plausible that future simulations of the full 
FTGPE (||) or approximations to it could be used to 
quantitatively measure the shift in the critical temper- 
ature for the homogeneous Bose gas when the lowest 
energy modes are sufficiently classical. However, the 
terms coupling the FTGPE to the effective bath ?7(x) 
may be difficult to implement computationally, and at 
the present time we are unsure how to proceed in this 
direction. 



not have any results to compare with for the strongly 
interacting regime where the non-perturbative method 
was used. We can only conclude that the temperatures 
extracted using this method agree numerically with the 
other two methods in the weakly interacting regime, and 
that the values obtained seem reasonable and basis inde- 
pendent elsewhere. We intend to test this method further 
in the future using a numerical "ideal gas thermometer" . 

In this section we move on to consider how other sys- 
tem properties such as condensate fraction, specific heat, 
and vorticity vary with the temperature T. 



A. Condensate fraction 

It is usual when considering how the condensate frac- 
tion varies with the other properties of the system to plot 
it against temperature, rather than against energy as we 
have done in Fig. ^. We are now in a position to present 
this data, and it is displayed in Fig. ^ We can see that a 
major effect of increasing the nonlinearity is to increase 
the condensate fraction at any given temperature. This 
can be understood in the Bogoliubov regime by consid- 
ering the shape of the dispersion relation. 

The Bogoliubov dispersion relation Eq. (^ ) shows that 
for a given condensate fraction, a larger value of Cni will 
result in an increase in the energy of any mode k relative 
to the condensate. This leads directly to the observation 
that for a fixed condensate fraction, an increase in the 
nonlinearity must lead to an increase in the temperature. 
However, as {Nq)/N ^ in the transition region, the 
energy-momentum relationship tends towards the ideal 



B. Specific heat 

In Fig. |l^(a) we plot the energy of the simulations due 
to excited states (E — Eq) versus temperature, where 
Eq = Cni/2 is the energy of the system at T = 0. We can 
see that at low temperatures for all interaction strengths 
this is a straight line, with a slope of about 13996 — 
the number of excited modes in the system. This is as 
expected — the average energy contained in a given mode 
k is 

(E,) EE = « T, (43) 

when there is a condensate present and p, —^ 0^ . At 
higher temperatures, however, the energy rises above the 
equipartition prediction for non-zero interaction strength 
Cni, and this is an indication that there are no longer 
independent modes in the system. 

The derivative of this curve with respect to tempera- 
ture gives the specific heat, and this quantity is plotted 
in Fig. |l^(b). The messy nature of this plot is due to 
small uncertainties in the measured temperature which 
are amplified when the temperature difference between 
successive simulations is calculated. However, the plot 
does display an interesting feature. For non-zero interac- 
tion strength, the specific heat appears to reach a peak at 
the transition temperature, and the height of this peak 
increases with the value of Cni — somewhat reminiscent of 
the lambda transition in superfluid helium. Once again, 
further simulations and more accurate determination of 
the temperature are required for quantitative investiga- 
tion of this effect. This will be the subject of future work. 
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FIG. 10: Graphs relating to the specific heat of the Bose gas 
in the PGPE model. The open circles are for Cni = 10000, 
crosses for Cni ~ 2000, solid dots for Cni = 500, and the 
solid line is for the ideal gas. (a) Plot of the energy versus 
temperature for all four interaction strengths considered, (b) 
Plot of the specific heat for all four interaction strengths. 



The role of vortices 



A further quantity of interest is the vorticity of the 
system in equilibrium. It has been argued that vortices 
may be important in the superfluid transition of *He, re- 
ducing the superfluid density near the transition point 
p5[ . With this in mind, we have studied the presence 
of vortex lines in our simulations. Recently Berloff and 
Svistunov ||2^ have considered the evolution of topologi- 
cal defects in the evolution of a Bose gas from a strongly 
non-equilibrium state. 

A vortex is a topological excitation, characterised in a 



(j) VArg[i/'(x)] -^1 = 27 



(44) 



where C is a closed contour, and n is a non-zero integer, 
the sign of which indicates the circulation of the vortex. 
The continuous variation of the phase from zero to 2m: 
around such a contour implies that there must be a dis- 
continuity in the phase within the loop. The only way 
that this can be physical is for the wave function to have 
zero amplitude at the spatial position of the phase sin- 
gularity. 

In a two-dimensional wave function the centre of vor- 
tices are zero-dimensional points, and they can be easily 
counted to give a measure of the vorticity of the sys- 
tem. However, in three dimensions vortices form lines 
and rings, and the equivalent quantity of the 2D mea- 
sure of vorticity would be to calculate the length of all 
vortex structures in the wave function. This would be a 
somewhat complicated procedure numerically, and so we 
have devised a different technique. 

We increase the spatial resolution of our wave functions 
to be 128 X 128 x 128 points, so that the grid spacing is 
smaller than the vortex healing length ^, defined by 



2mt 



= noUo 



(45) 



We do this by extending the wave function in fc-space, and 
then Fourier transforming to real space. This does not 
require any extra information, as for k > 15 x 2'k/L we 
have Ck = 0. We then count the number of vortex lines 
passing through every xy plane, and take the average over 
all planes. It seems that this is a reasonable measure of 
the vorticity of the wave function, and it should be similar 
to the measurement of the length of the vortex structures 
discussed above. 

We have analysed the data from the simulations using 
this procedure. We find that when the energy of the sim- 
ulation is sufficiently high that there are vortices present, 
the time evolution of the vorticity is a good indicator for 
when the system reaches equilibrium. As is the case for 
the condensate population, the vorticity tends to an equi- 
librium value which fluctuates by a small amount (much 
smaller than the fluctuations in the condensate popula- 
tion). 

A plot of the vorticity against system energy is shown 
in Fig. |l^(a) for the Cni = 10000 simulation (the curves 
are qualitatively similar for the other nonlinearities). We 
see that there is a minimum energy required for vortices 
to be present in the system at equilibrium. Also, as we 
reach this energy the plot of condensate occupation ver- 
sus energy appears to dip. This same behaviour is ob- 
served for the Cni = 500 and 2000 cases, but it occurs at 
a higher condensate fraction, and is not as pronounced. 
There is no corresponding departure from linearity in the 
ideal gas case, as was seen in Fig. ^(b). 
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FIG. 11: The presence of vortices in the simulations, (a) 
A plot of vorticity for the Cni = 10000 simulation series. 
The number of vortex lines per plane are indicated by open 
circles with the scale on the left vertical axis, and the con- 
densate fraction by dots with the scale on the right vertical 
axis, (b) The number of vortex lines per plane plotted against 
temperature for all three simulation series. Open circles are 
Cni = 10000, crosses are Cni = 2000, and dots are Cni = 500. 



A plot of the number of vortex lines versus temperature 
for all the simulations is shown in Fig. pT|(b), and this dis- 
plays a large increase in the vorticity near the transition 
temperature for Cni = 10000. Even the Cni = 2000 case 
appears to show a small jump in this region. A more 
in-depth analysis of this behaviour will be carried out in 
a subsequent extension of this work. 

Finally, a three-dimensional visualisation of the net- 
work of vortex lines is shown in Fig. |l^ for three simula- 
tion energies for the Cni = 10000 simulations. Each point 
corresponds to where a vortex line was detected in the 
horizontal planes, and for the lowest two energies several 
vortex rings are clearly visible. 
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FIG. 12: A visualisation of the "vortex tangle" in equilibrium 
for the case of Cni = 10000. (a) E = 7000, (b) E = 8000, 
(c) E = 9000. Each point corresponds to where a vortex line 
was detected in the horizontal plane. Several vortex rings are 
visible in the figures. 



16 



VIII. CONCLUSIONS 

We have presented what we believe is compelhng ev- 
idence that the projected Gross-Pitaevskii equation is 
a good approximation to the dynamics of the classical 
modes of a Bose gas. We have described how to carry out 
the projection technique in the homogeneous case with 
periodic boundary conditions, and have shown that start- 
ing with a randomised wave function with a given en- 
ergy, the projected GPE evolves towards an equilibrium 
state. We have analysed the numerical data in terms of 
quadratic Bogoliubov theorj^and also the gapless, finite 
temperature theory of Ref. Q in the classical limit. We 
have found that both the occupation and energies of the 
quasiparticles agree quantitatively with the predictions 
when these theories are valid. 

Outside the range of perturbation theory we have pro- 
posed another technique that has allowed us to determine 



a temperature for the PGPE simulations in equilibrium. 
This method agrees with the perturbative methods when 
they are valid. Using this definition, we have found that 
increasing the nonlinearity Cni leads to both an increase 
in the transition temperature, and in the specific heat of 
the system at the critical point. We have also presented 
evidence that suggests vortices may play some role in 
the transition. The projected GPE is a simple equation 
but it appears to describe very rich physics, only some of 
which we have considered here. 
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